function cholesky_applications()
% CHOLESKY_APPLICATIONS Cholesky分解的应用
% 
% 展示Cholesky分解在各种问题中的应用

fprintf('=== Cholesky分解应用 ===\n\n');

%% 1. 线性系统求解
fprintf('1. 线性系统求解\n');
fprintf('目标：使用Cholesky分解求解正定系统\n\n');

% 构造正定矩阵
A1 = [4, 2, 1; 2, 5, 3; 1, 3, 6];
b1 = [1; 2; 3];

fprintf('正定矩阵 A1 =\n');
disp(A1);
fprintf('右端向量 b1 = [%g, %g, %g]''\n', b1);

% 验证正定性
eigenvals1 = eig(A1);
fprintf('特征值: [%.4f, %.4f, %.4f]\n', eigenvals1);
fprintf('正定性: %s\n', mat2str(all(eigenvals1 > 0)));

% Cholesky分解
L1 = chol(A1, 'lower');
fprintf('Cholesky因子 L1 =\n');
disp(L1);

% 验证分解
A1_reconstructed = L1 * L1';
fprintf('重构误差: %.2e\n', norm(A1 - A1_reconstructed, 'fro'));

% 求解线性系统 A1*x = b1
% 步骤1: 求解 L1*y = b1
y1 = L1 \ b1;
% 步骤2: 求解 L1'*x = y1
x1 = L1' \ y1;

fprintf('解向量 x1 = [%.6f, %.6f, %.6f]''\n', x1);

% 验证解
residual1 = norm(A1 * x1 - b1);
fprintf('残差: %.2e\n', residual1);

% 与直接求解比较
x1_direct = A1 \ b1;
fprintf('直接求解结果: [%.6f, %.6f, %.6f]''\n', x1_direct);
fprintf('解的差异: %.2e\n', norm(x1 - x1_direct));

%% 2. 最小二乘问题的正规方程
fprintf('\n2. 最小二乘问题的正规方程\n');
fprintf('目标：使用Cholesky分解求解正规方程\n\n');

% 过定系统
m = 10; n = 4;
A2 = randn(m, n);
x_true = [1; -2; 3; 0.5];
noise = 0.1 * randn(m, 1);
b2 = A2 * x_true + noise;

fprintf('过定系统: %d×%d 矩阵\n', m, n);
fprintf('真实解: [%.2f, %.2f, %.2f, %.2f]''\n', x_true);

% 正规方程: A'A*x = A'b
AtA = A2' * A2;
Atb = A2' * b2;

fprintf('正规方程矩阵 A''A 的条件数: %.2e\n', cond(AtA));

% 检查正定性
eigenvals_AtA = eig(AtA);
fprintf('A''A 的最小特征值: %.2e\n', min(eigenvals_AtA));

if min(eigenvals_AtA) > 1e-12
    % Cholesky分解求解
    L2 = chol(AtA, 'lower');
    y2 = L2 \ Atb;
    x2_chol = L2' \ y2;
    
    fprintf('Cholesky求解结果: [%.4f, %.4f, %.4f, %.4f]''\n', x2_chol);
    
    % 与QR分解比较
    [Q2, R2] = qr(A2, 0);
    x2_qr = R2 \ (Q2' * b2);
    
    fprintf('QR分解求解结果: [%.4f, %.4f, %.4f, %.4f]''\n', x2_qr);
    fprintf('两种方法的差异: %.2e\n', norm(x2_chol - x2_qr));
    
    % 残差比较
    residual_chol = norm(A2 * x2_chol - b2);
    residual_qr = norm(A2 * x2_qr - b2);
    fprintf('残差比较: Cholesky = %.2e, QR = %.2e\n', residual_chol, residual_qr);
else
    fprintf('A''A 不是正定的，无法使用Cholesky分解\n');
end

%% 3. 协方差矩阵的处理
fprintf('\n3. 协方差矩阵的处理\n');
fprintf('目标：分析和处理协方差矩阵\n\n');

% 生成相关数据
n_samples = 100;
n_vars = 4;
rng(42);

% 构造具有特定相关结构的数据
L_true = [2, 0, 0, 0; 1, 1.5, 0, 0; 0.5, 0.8, 1, 0; 0.2, 0.3, 0.5, 0.8];
X3 = randn(n_samples, n_vars) * L_true';

fprintf('数据维度: %d×%d\n', size(X3));

% 计算样本协方差矩阵
X3_centered = X3 - mean(X3);
Cov3 = (X3_centered' * X3_centered) / (n_samples - 1);

fprintf('样本协方差矩阵 =\n');
disp(Cov3);

% 检查正定性
eigenvals_cov = eig(Cov3);
fprintf('协方差矩阵特征值: [%.4f, %.4f, %.4f, %.4f]\n', eigenvals_cov);
fprintf('最小特征值: %.4f\n', min(eigenvals_cov));

if min(eigenvals_cov) > 1e-10
    % Cholesky分解
    L3 = chol(Cov3, 'lower');
    fprintf('协方差矩阵的Cholesky因子 L3 =\n');
    disp(L3);
    
    % 生成具有相同协方差的新数据
    Z3 = randn(50, n_vars);  % 标准正态数据
    X3_new = Z3 * L3' + mean(X3);  % 变换到目标分布
    
    % 验证新数据的协方差
    X3_new_centered = X3_new - mean(X3_new);
    Cov3_new = (X3_new_centered' * X3_new_centered) / (size(X3_new, 1) - 1);
    
    fprintf('新数据的协方差矩阵 =\n');
    disp(Cov3_new);
    fprintf('协方差差异: %.2e\n', norm(Cov3 - Cov3_new, 'fro'));
else
    fprintf('协方差矩阵不是正定的\n');
end

%% 4. 多元正态分布的模拟
fprintf('\n4. 多元正态分布的模拟\n');
fprintf('目标：使用Cholesky分解生成多元正态随机数\n\n');

% 目标分布参数
mu4 = [1; 2; -1];
Sigma4 = [2, 0.5, 0.2; 0.5, 1.5, 0.8; 0.2, 0.8, 1];

fprintf('目标均值: [%.1f, %.1f, %.1f]''\n', mu4);
fprintf('目标协方差矩阵 =\n');
disp(Sigma4);

% 检查正定性
eigenvals_sigma = eig(Sigma4);
fprintf('协方差矩阵特征值: [%.4f, %.4f, %.4f]\n', eigenvals_sigma);

if all(eigenvals_sigma > 0)
    % Cholesky分解
    L4 = chol(Sigma4, 'lower');
    fprintf('Cholesky因子 L4 =\n');
    disp(L4);
    
    % 生成多元正态随机数
    n_samples4 = 1000;
    Z4 = randn(n_samples4, 3);  % 标准正态
    X4 = Z4 * L4' + mu4';       % 变换到目标分布
    
    % 验证生成的数据
    mu4_sample = mean(X4)';
    Sigma4_sample = cov(X4);
    
    fprintf('样本均值: [%.4f, %.4f, %.4f]''\n', mu4_sample);
    fprintf('样本协方差矩阵 =\n');
    disp(Sigma4_sample);
    
    fprintf('均值误差: %.2e\n', norm(mu4 - mu4_sample));
    fprintf('协方差误差: %.2e\n', norm(Sigma4 - Sigma4_sample, 'fro'));
    
    % 概率密度函数计算
    test_point = [1.5; 2.2; -0.8];
    diff = test_point - mu4;
    
    % 使用Cholesky因子计算
    y = L4 \ diff;
    log_det_sigma = 2 * sum(log(diag(L4)));  % log(det(Sigma)) = 2*sum(log(diag(L)))
    log_pdf = -0.5 * (3*log(2*pi) + log_det_sigma + y'*y);
    
    fprintf('测试点 [%.1f, %.1f, %.1f] 的对数概率密度: %.4f\n', ...
            test_point, log_pdf);
else
    fprintf('协方差矩阵不是正定的\n');
end

%% 5. 优化问题中的应用
fprintf('\n5. 优化问题中的应用\n');
fprintf('目标：求解二次规划问题\n\n');

% 二次规划问题: min (1/2)x'Hx + c'x
% 其中 H 是正定矩阵
H5 = [2, 1; 1, 3];
c5 = [-1; -2];

fprintf('Hessian矩阵 H5 =\n');
disp(H5);
fprintf('线性项系数 c5 = [%.1f, %.1f]''\n', c5);

% 检查正定性
eigenvals_H = eig(H5);
fprintf('Hessian特征值: [%.4f, %.4f]\n', eigenvals_H);

if all(eigenvals_H > 0)
    % 最优解: x* = -H^(-1)*c
    % 使用Cholesky分解求解
    L5 = chol(H5, 'lower');
    y5 = L5 \ (-c5);
    x5_opt = L5' \ y5;
    
    fprintf('最优解: x* = [%.6f, %.6f]''\n', x5_opt);
    
    % 验证：梯度应该为零
    gradient = H5 * x5_opt + c5;
    fprintf('最优点处的梯度: [%.2e, %.2e]''\n', gradient);
    
    % 计算最优值
    opt_value = 0.5 * x5_opt' * H5 * x5_opt + c5' * x5_opt;
    fprintf('最优目标函数值: %.6f\n', opt_value);
    
    % 与直接求解比较
    x5_direct = -H5 \ c5;
    fprintf('直接求解结果: [%.6f, %.6f]''\n', x5_direct);
    fprintf('解的差异: %.2e\n', norm(x5_opt - x5_direct));
else
    fprintf('Hessian矩阵不是正定的\n');
end

%% 6. 数值稳定性分析
fprintf('\n6. 数值稳定性分析\n');
fprintf('目标：分析Cholesky分解的数值稳定性\n\n');

% 测试不同条件数的正定矩阵
condition_numbers = [1e2, 1e6, 1e10, 1e12];
fprintf('条件数        分解成功    重构误差      求解误差\n');
fprintf('--------      --------    --------      --------\n');

for cond_target = condition_numbers
    % 构造指定条件数的正定矩阵
    n = 4;
    eigenvals_target = [1, 1/sqrt(cond_target), 1/cond_target^(1/3), 1/cond_target];
    V_rand = orth(randn(n, n));
    A_test = V_rand * diag(eigenvals_target) * V_rand';
    
    % 确保对称性
    A_test = (A_test + A_test') / 2;
    
    try
        % Cholesky分解
        L_test = chol(A_test, 'lower');
        decomp_success = true;
        
        % 重构误差
        A_reconstructed = L_test * L_test';
        reconstruction_error = norm(A_test - A_reconstructed, 'fro') / norm(A_test, 'fro');
        
        % 求解误差
        b_test = randn(n, 1);
        x_chol = L_test' \ (L_test \ b_test);
        x_direct = A_test \ b_test;
        solve_error = norm(x_chol - x_direct) / norm(x_direct);
        
        fprintf('%.0e      %8s    %.2e      %.2e\n', ...
                cond_target, 'Yes', reconstruction_error, solve_error);
        
    catch
        fprintf('%.0e      %8s    %8s      %8s\n', ...
                cond_target, 'No', 'N/A', 'N/A');
    end
end

%% 7. 修正Cholesky分解
fprintf('\n7. 修正Cholesky分解\n');
fprintf('目标：处理非正定矩阵\n\n');

% 构造非正定矩阵
A7 = [2, 1, 1; 1, 1, 1; 1, 1, 1];  % 奇异矩阵
fprintf('非正定矩阵 A7 =\n');
disp(A7);

eigenvals7 = eig(A7);
fprintf('特征值: [%.4f, %.4f, %.4f]\n', eigenvals7);
fprintf('最小特征值: %.4f (≤ 0，非正定)\n', min(eigenvals7));

% 修正方法1：添加对角扰动
tau = 1e-3;
A7_modified1 = A7 + tau * eye(size(A7));
fprintf('方法1：添加对角扰动 τ = %.0e\n', tau);

try
    L7_mod1 = chol(A7_modified1, 'lower');
    fprintf('修正后Cholesky分解成功\n');
    
    % 验证修正后的正定性
    eigenvals7_mod1 = eig(A7_modified1);
    fprintf('修正后最小特征值: %.4f\n', min(eigenvals7_mod1));
catch
    fprintf('修正后仍无法进行Cholesky分解\n');
end

% 修正方法2：基于特征值分解的修正
[V7, D7] = eig(A7);
eigenvals7_corrected = max(eigenvals7, tau);
A7_modified2 = V7 * diag(eigenvals7_corrected) * V7';

fprintf('方法2：特征值修正\n');
fprintf('修正后特征值: [%.4f, %.4f, %.4f]\n', eigenvals7_corrected);

try
    L7_mod2 = chol(A7_modified2, 'lower');
    fprintf('特征值修正后Cholesky分解成功\n');
catch
    fprintf('特征值修正后仍无法进行Cholesky分解\n');
end

%% 总结
fprintf('\n=== Cholesky分解应用总结 ===\n');
fprintf('1. 线性系统求解：效率高，数值稳定\n');
fprintf('2. 最小二乘正规方程：需要注意条件数\n');
fprintf('3. 协方差矩阵处理：自然的正定结构\n');
fprintf('4. 多元正态分布：高效的随机数生成\n');
fprintf('5. 优化问题：二次规划的标准方法\n');
fprintf('6. 数值稳定性：取决于矩阵条件数\n');
fprintf('7. 修正技术：处理非正定情况\n\n');

fprintf('Cholesky分解的优势：\n');
fprintf('• 计算量约为LU分解的一半\n');
fprintf('• 数值稳定性好\n');
fprintf('• 保持正定性\n');
fprintf('• 存储需求少\n\n');

fprintf('使用注意事项：\n');
fprintf('• 仅适用于正定矩阵\n');
fprintf('• 高条件数时精度下降\n');
fprintf('• 需要检查正定性\n');
fprintf('• 必要时进行修正\n\n');

end